Semiclassical Theory of Integrable and Rough Andreev Billiards 
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We study the effect on the density of states in mesoscopic ballistic billiards to which a super- 
conducting lead is attached. The expression for the density of states is derived in the semiclassical 
S-matrix formalism shedding insight into the origin of the differences between the semiclassical the- 
ory and the corresponding result derived from random matrix models. Applications to a square 
billiard geometry and billiards with boundary roughness are discussed. The saturation of the quasi- 
particle excitation spectrum is related to the classical dynamics of the billiard. The influence of 
weak magnetic fields on the proximity effect in rough Andreev billiards is discussed and an analyti- 
cal formula is derived. The semiclassical theory provides an interpretation for the suppression of the 
proximity effect in the presence of magnetic fields as a coherence effect of time reversed trajectories, 
similar to the weak localisation correction of the magneto-resistance in chaotic mesoscopic systems. 
q , The semiclassical theory is shown to be in good agreement with quantum mechanical calculations. 
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^ \ PACS numbers: 05.45+b, 74.50+r, 74.80Fp 
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A superconductor in proximity to a normal conductor affects the spectral density of quasiparticle excitations in the 
conductor. Recent technological advances in building very clean conductors of mesoscopic size have led to consider 
this proximity effect not only in the dirty disorder limit, where it is known to play an important role in many transport 
properties but also in ballistic mesoscopic samples J5|-|| in proximity to a bulk superconductor. These systems 

are often modelled by ballistic billiards. Billiard shaped structures connected to a superconductor have been coined 
Andreev billiards H. It was shown in [^,^) that the form of the quasiparticle excitation spectrum in the Andreev 
billiard depends crucially on whether its classical dynamics is integrable or chaotic. A semiclassical interpretation of 
these results based on the Eilenberger Green's function jl(J was given in [Q. The semiclassical approach has also been 
. extended to mesoscopic samples with a mixed classical phase space [gj . 

The aim of this contribution is twofold: We first present a derivation of the semiclassical result for the quasiparticle 
excitation spectrum based on the semiclassical scattering matrix approach as pioneered by Smilansky and co-workers 
|T[| , |r2[ which allows for a transparent physical interpretation of the result of j7j] in terms of multiple Andreev scattering 
events [ p~3| . Secondly, we discuss a number of features of the quasiparticle excitation spectrum in an integrable square 
billiard and in a square billiard with surface roughness. These features include the saturation of the quasiparticle 
^ . excitation spectrum in the square billiard and a semiclassical explanation of the effect a weak magnetic field has on 
03 the quasiparticle excitation spectrum in billiards with surface roughness. 

The semiclassical expression for the quasiparticle density of states will be derived in Section 2. The strength of this 
approach lies in the fact that it allows to see on which level approximations enter scmiclassically. The semiclassical 
theory of chaotic Andreev billiards predicts an exponential suppression of the density of states near the Fermi energy 
[0,^) . In contrast a random matrix modelling leads to the result of a gap in the density of states in an energy 
interval above the Fermi energy. The scattering matrix approach elucidates one possible origin of the differing results. 

In Sec. 3.1 we discuss the quasiparticle excitation spectrum of a square Andreev billiard in detail. While previous 
papers concentrated on the linear rise of the spectrum above the Fermi energy l^-Q we focus on the saturation of the 
excitation spectrum at higher energies. We show that the saturation can be related to the probability distribution of 
short classical paths hitting the superconducting parts of the billiard boundaries. In Sec. 3.2 we discuss results for 
square billiards with additional surface roughness. Surface roughness is modelled as isotropic scattering of electrons 
and holes from the normal parts of the billiard boundary. An exponential suppression of the density of states near the 
Fermi energy is observed in complete analogy to chaotic billiards ||. Finally in Sec. 3.3 the effect of an additional weak 
magnetic field on the quasiparticle excitation spectrum in the rough billiard is considered. We predict an enhanced 
density of states compared to the field free case for energies near Ef and give analytical semiclassical expressions. 
The semiclassical theory provides a clear interpretation of this result as an effect of destructive interference between 
reversed paths in the presence of a magnetic field. Again the semiclassical predictions are in good agreement with 
quantum mechanical calculations. 
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2 Theory 



In a quasi-classical picture the effect of a superconducting lead coupled to a mesoscopic billiard manifests itself 
in the process known as Andreev reflection Jl3| ] : At the superconducting parts of the billiard boundary electron-like 
quasi-particles are retro-reflected with opposite velocities as holes and vice versa (see Fig. [l], and Refs. HH for a 
detailed description) . In contrast electrons and holes are specularly reflected at those parts of the billiard boundary 
which are not in contact to the superconductor if the billiard boundary is smooth. For a billiard with surface roughness 
we will take the possibility of isotropic scattering off the normal walls of the billiard into account. This is also indicated 
in Fig. |. 

The local quasiparticle density of states (DOS) d(E, r) of an Andreev billiard is defined as the density of states at 
positive energy E above the Fermi energy Ep = 0, weighted with the corresponding electron-like component of the 
wave function at point r. It is given by |?J 

A denotes the area of the billiard and vf denotes the Fermi velocity. It is assumed that E is much smaller than the pair 
potential A(r) in the superconducting lead and that the lead supports a large number N ^> 1 of classically allowed 
transverse channels. In the billiard the pair potential A(r), which couples electron and hole like states, vanishes 
identically. Due to the condition E <C A electrons in the billiard which hit the superconducting lead get reflected 
as holes thus quantum mechanically forming electron-hole quasiparticle states. L(cj>) is the length of the trajectory 
which passes the point r in the billiard under an angle (f> between two successive bounces with the superconducting 
boundary and djv = mA/ {2irfi 2 ) is the average density of states in the isolated billiard. We further assume a perfectly 
transmitting boundary between the billiard and the lead if the lead is in the normal state (no potential difference), 
which has the effect that the probability for Andreev reflection equals one if the lead is in the superconducting state. 
(A situation including the probability for Andreev reflection as well as normal reflection between SN boundaries was 
taken into account in a somewhat different context and geometry e.g. in |]l5,|l6| ) 



The total quasiparticle DOS in the billiard is obtained by integrating (l|) over the area of the billiard system. 
The integral over the billiard area can be converted into an integral over initial starting angles a and positions y of 
trajectories along the lead of width w. The resulting expression is 

io 1 L(y,a) ^ 

d(E) = d -fjdyJ d(sina) J ds^J (^g^ - (»+ ^r) , (2) 



o-i o 



where s is the local variable measuring the length along a trajectory. If one additionally assumes an ergodic distribution 
of trajectories in the initial conditions y and sin a on the boundary the expression for the density of states can be 
rewritten in terms of the probability distribution P(L) of trajectories of length L between two successive bounces 
with the SN boundary. It is then 



(n + -)nj . (3) 

The assumption of an ergodic distribution is justified in the case of rough billiards when the time t CTg on which the 
classical motion in the billiard becomes ergodic is much smaller than the mean escape time r csc of a particle from 
the billiard into the SN-lead. It is also a good approximation for the integrable square billiard with a large number 
of open channels N ^S> 1 because then the distribution of initial conditions in sin(? is quasi-continuous. For a small 
number of channels it is more appropriate to resort to the continued fraction evaluation of (^J) described in Sec. Ilia, 
which takes the quantization of angles explicitely into account. 

We now aim to derive (|^) from the semiclassical scattering matrix approach for quantisation ]l2[ ] and to understand 
the approximations which lead to (^|) on the semiclassical level. The scattering matrix approach for quantisation of 
a billiard system which is opened via a lead starts with a trace formula for the density of states d{E) which involves 
two contributions: The first term dn(E) is the resonance density of states in the corresponding billiard opened via 
the lead. The second contribution takes the coupling to the superconductor into account and is expressed as a sum 
over traces of powers of the scattering matrix which relates incoming and outgoing transverse modes in the lead. The 
result is: 
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d(E) = d R (E) - - lim Im— lndet [1 - S A S N (E + if])} 

7T rj^O oE 



(4) 



where the total scattering matrix for the Andreev billiard is composed of a product of a 2N times 2N normal scattering 
matrix Sn(E) and an Andreev scattering matrix Sa of equal dimension. The two scattering matrices are given as 



S N (E) 



o 

*(-E) 



and 



S A = 







(•5) 



The normal scattering matrix S n has block diagonal structure where the two N times N blocks describe electron and 
hole scattering between the channel modes. The Andreev scattering matrix Sa couples electrons and holes at the 
SN-interface. The sub-diagonal elements are N times N unit matrices with an additional phase of (—i). We neglect 
the weak energy dependence of the Andreev scattering matrix which is valid in the deep sub-gap regime E <C A 
flilpol . In terms of the electron and hole scattering matrix the density of states takes the form 



1 °° ( 1 \ m f) 
d{E) = dji(E) - -Im £ i 1 ^Tr— [S(E)S*(-E)] m , 



(6) 



where (||) has additionally been expanded into a sum over traces of powers of the scattering matrix. An averaged 
density can be obtained by different procedures: One can average d(E) over a classically small interval of Fermi 
energies Ep or one can average d(E) over different realizations of the rough billiard (with fixed area A). In any case 
averaging of the resonance density gives the average density djy of the isolated billiard, and the average quasiparticle 
density is 



d av {E) = d N (E) 



1 



-Im 



E 



-ly 



■Tr-^([S(E)S*(-E)] m ) 



(7) 



Averaging is denoted by brackets (•••). Equation (Q) is the starting point for the derivation of expression (^) within 
the semiclassical scattering matrix approach. 

Scmiclassically the elements S nn > (E) of the electron scattering matrix are expressed as a sum over classical orbits. 
Each orbit contributes with an amplitude Aj and a phase Sj. The resulting expression is [pT| , ^9[ 



S nn > (E) = Aj (n — >■ n') exp 



1 jr 



(8) 



where Uj is an additional integer Maslov index. The amplitude pre-factor is given explicitely by 



Aj(n^n') = — 



2ir 



^1/2 


di'(E) 




de 



-1/2 



(9) 



where I' — Tin' is the action of the final transverse motion in the lead. Only those paths contribute which enter the 
billiard at the SN-boundary with fixed quantised angle ±sin# = mr / '(kpW) and return to the boundary with angle 
±sin#' = n 'it / '(kpW). In terms of the initial position y along the SN boundary and the angle 8' with which the 
trajectory returns the amplitude can be written as 



Aj{n 



TT 
2kF~ 



dy 



d(sm0') 



1/2 



(10) 



To proceed further we resort to the physical picture of Andreev reflection of electrons into holes and vice versa at the 
SN-boundary. We first observe that the traces contain products of alternating electron scattering matrices S(E) and 
hole scattering matrices S*(—E). The energy E above the Fermi level is large with respect to the mean level density 
8 = d^ 1 of the isolated billiard, but is classically small. Similar to the semiclassical evaluation of density-density 
correlator ]26|-|28[] we expand the phase around the Fermi energy as Sj (±E) ~ Sj (0) ± ETj (0) where Tj is the return 
time of the orbit to the SN-interface. Additionally the amplitudes are only slowly varying functions of the energy and 
are evaluated at the Fermi energy. 
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We demonstrate how to evaluate the sum over traces of products of the scattering matrices for the n = 1 term. 
Generalisation to higher order terms is straightforward. Matrix elements of products of an electron and a hole 
scattering matrix have the form 



[S(E)S*(-E)] nn , = J2J2 A o( n -» n ") A U n " -» « ex P 

n" j,k 



exp 



(11) 



Upon averaging [S(E)S*(— E)] nn i over the Fermi energy or different realizations of boundary roughness in the case 
of rough billiards only diagonal terms j = k contribute to the sum. The diagonal approximation is further justified 
by the physical picture of Andreev reflection, which means that the reflected hole orbit retraces the electron orbit 
and vice versa. However it must be emphasised that Andreev reflection is exactly fulfilled only at the Fermi energy 
Ep = 0, and in an exact treatment deviations from Andreev reflection at finite energies have be taken into account. 

The diagonal approximation implicates that the initial and final channel indices are equal, n = n' , and the product 
matrix (|ll|) is diagonal: 

[S(E)S*(-E)] nn , = S nnr J2J2\ A ^ n ^ n ")\ 2 ex ^ ( 2 i ET ^j ■ ( 12 ) 

n" j \ 1 J 

In the semiclassical limit the summation over intermediate channels n" can be transformed into an integral over 
angles: ^fw/tt d(sm9'). Using the expression ( |Io| ) for the amplitudes one arrives at the final expression 

1 r w i 

[S(E)S*(-E)] nn> =6 nn ,— / dyexp(2-ET(y)). (13) 
2w J h 

Taking the trace amounts to another integration over the angle in the semiclassical approximation and one has 

Tt[S(E)S*(-E)} = ^-J d(wx6) J dyexp(2-ET(y)). (14) 

Within the diagonal approximation traces of higher powers of products of electron and hole scattering matrices are 
easily shown to be 

r Tr[S{E)S*{-E)] m = 1 f d( S in6) dyexp(2^mET(y)) . (15) 
Using the trace formula (0) this gives 

d av (E) =dAT + -%- VM) m / dy I d(sm6)Tcos(^E). (Hi) 



2n 2 h 



^(-ir/ dy d(sind)Tco S (^-E). 



If now Poissonian summation is used the equivalence of ( |l6| ) and the Bohr-Sommerfeld like expression (g) for the 
average density of states in the Andreev billiard, which was derived on the basis of the Eilenberger equation for the 
Green's function ||, can easily be seen. 

The scattering matrix approach gives a clear and intuitive interpretation of how the coupling of the billiard to 
superconducting leads modifies the average density of states: It is given as a sum of the average density of states 
of the isolated billiard (the Weyl term) plus a sum of multiple Andreev reflections of electron into hole states. Off- 
diagonal corrections arise due to the fact that the condition of Andreev reflection is exactly fulfilled only at the Fermi 
energy. The presence of non-diagonal contributions points towards a possible explanation for the difference between 
the random matrix result for a chaotic Andreev billiard J^,^|, which predicts a gap in the quasiparticle excitation 
spectrum for billiards with a chaotic classical phase space, and the semiclassical theory which leads to an exponential 
suppression |,|. Pairs of non-identical trajectories k ^ j which follow the same path along a segment in real space 
before they depart due to slightly different initial conditions and the chaotic dynamics may have similar actions and 



therefore also survive the averaging of products of electron and hole scattering matrices, Eq. (11) and powers of it pi| . 
An extended diagonal approximation which takes the contribution of paths into account which are not exactly related 
by time-reversal symmetry has recently been proposed in the context of weak localisation |^,^3|. The influence of 
off-diagonal contributions arising in the semiclassical approach to the density of states in Andreev billiards remains a 
topic to be investigated. 
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Finally it should be emphasised that averaging has been performed on classically small scales. Fluctuations of 
d av (E) can still appear on classical energy scales due to fluctuations in the length distribution probability P(L). In 
the following we will skip the index. d(E) then also denotes the density of states after the above described averaging. 



3 Results 

In the following we present results of quantum mechanical calculations for the average density of states in Andreev 
billiards and compare them with the semiclassical theory. We will focus on the difference between results for a square 
billiard which is integrable and a square billiard with rough boundaries. As an effect of rough boundaries electron and 
hole trajectories can be scattered into arbitrary directions when they hit the normal boundary. We will then consider 
the effect of a weak magnetic field on the density of states and derive a semi-analytical expression for the magnetic 
flux and energy dependence of the density of states. 

To numerically model the structure of Fig.[j] we consider a ballistic normal region connected to a clean supercon- 
ductor. For our numerical calculation we use a tight-binding version of the Bogoliubov-de Gennes Hamiltonian: 

H, A \(u\ _ E (u\ 



A* -H J \ v J \v 

In this equation Ho = J2i ~ I*) 01 ^ s the standard single-particle Anderson model with (ij) denoting 

pairs of nearest neighbour sites and A = J^. |i)Aj(i| is the superconducting order parameter. The billiard has width 
M and length L (in units of the lattice constant a). The coupling to the superconductor is of width w. With the 
exception of the billiard boundary, the diagonal matrix elements e; = eo, with eo chosen to ensure that the Fermi 
level is away from the van Hove singularity in the band centre (eo = 0). To model the boundary roughness, for sites 
on the boundary we randomly choose = 10 4 or = eo with equal probability. In the absence of a magnetic field, 
the off-diagonal matrix elements t, which determine the width of the energy band (the band- width is 8t), are set to 
1 throughout the system. The effects of a magnetic field are modelled by including a Peierls' phase factor into these 
off-diagonal elements in the billiard region. To compute the density of states, a numerical decimation technique is 
employed |Q . This has been used extensively over recent years to discuss transport and density of states properties 
of hybrid systems |25| . 

3.1 Square billiards 

The asymptotic length distribution of trajectories in square billiards typically follows a power law of the form 
P[L) ~ L~ 3 . This result was used in || to show that the density of quasiparticle states is a linear function of the 
energy of quasiparticle excitations for E <C A. The slope of the linear behaviour is related to the proportionality 
constant of the asymptotic length distribution. The above distribution however leads to a linear growth of the density 
of states at all energies and cannot reproduce the large energy limit d(E) —> g?at for the quasiparticle density. We 
show that the crossover to djv is related to a modified length distribution P(L) for small lengths. For widths w of the 
superconducting channel much smaller than the length a of the billiard the length distribution reaches approximately 
a plateau for small lengths. We introduce two different length distribution functions and compare numerical results 
to the analytical results derived from the smooth length distribution functions. 

It is useful to express the density of states in terms of scaled lengths and energies. Introducing the lengths 
Lt = ttA/w and Et — Tlvf/ '(2Lt) which are sometimes referred to as the Thouless length and the Thouless energy 
lengths I = L/Lt and energies e = E/Et are expressed in terms of these units. The meaning of Lt and Et 
will become clear for chaotic billiards (Sec. IIIB), where Lt is the average length of trajectories before they escape 
from the billiard. We use the same units of length and energy for the integrable square billiard in order to compare 
results with the rough billiard later. Expressed in scaled quantities the density of states has the form 

^- = ^1°° dW(l)l f>(j-(n+ i)^ (18) 



2' 



or equivalently 



i N e 2 ^ x 2 
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We introduce two different length distribution functions. The first length distribution function P C {L) approximates 
the probability by a constant for lengths I < l c : 



I < l c 
I > L 



The second distribution function has a smooth crossover to a constant probability at small I: 

Ps(l) ^ 



i 3 + i* 



(20) 



(21) 



We first discuss the length distribution function P c (l). For a given energy e there is a value no in the sum (|l9| ) 
given by n = el c /(2n) — 1/2 so that for n > n the algebraic tail of the length distribution for P c (l n ) applies while 
for n < Uq P c (l n ) is constant. The density of states is 



'JV 



"0 — 1 

73 E("+o) 



P 
c n=Q 



E 



(2tt) 3 ^ (n+i) 2 



which after summation and using the relation between no and e becomes 

2 



d(e) 
d N 



= C C 



2tt 2 /eZ c 1 

a? 1 2^ ~ 2 



2tt V2tt 



for e > 7r . 



(22) 



(23) 



■0' is the derivative of the Digamma function. This expression for the spectral density is valid for e > 7r, otherwise 
no = and the first part in the sum does not contribute. For e < it the result is 



d(e) _ n 
d N 4 



(24) 



which shows the known linear behaviour for small energies in square billiards Using *p'(x) 1/x a,s x — >oo the 
asymptotic behaviour for large energies is given by 



d(e) 



3C C 
2/ c 



as e 



oo . 



(25) 



The values of the free parameters are determined by the requirements that (a) the probability distribution must be 
normalised: J Q P{l)dl — 1 and (b) for large excitation energies the influence of the coupling to the superconductor 
vanishes and thus d(e) — > dN as e — > oo. This gives the values l c = 1 and C c = 2/3. It is seen that these requirements 
automatically determine the proportionality constant C c of the asymptotic power law tail of the length distribution 
( p0| ) and thus the slope C c 7r/4 of the linear behaviour ( p4| ) of the low energy spectral density. The value 7r/6 is 
somewhat smaller than the value 2/n given by the authors of Q who based their calculations on a numerically 
observed value of C = 8/ir 2 but did not take a crossover of P(l) to a flat distribution at some value l c into account. 
Note that one could adjust C independently to a given value if one introduced a three parameter function for the 
length distribution P(l), i.e. by giving a cut-off length below which P(?„) = 0. Such a cut-off length parameter 
would quite naturally be of order Z» « 2a. We will however not deal with this alternative here as when using a three 
parameter length distribution probability formulas become rather involved. 

The length distribution (20) is easy to handle and gives results which agree well with quantum mechanical calcula- 
tions. It leads however to an unphysical discontinuity of the derivative of d(e) at e — it. The second length distribution 
( f2~l|) we used is free of this feature. Proceeding in the same way as outlined above the spectral density of quasiparticle 
states is given by 



d N 



2tt ^ 



(26) 



The requirement of normalisation of P s (l) leads to C s = 3\/3l 2 / (2tt) . Further evaluation is only possible in the 
small and large energy limits. The small energy limit (el s <C 2tt) gives again the result (24) with C c replaced by 
C s . In the opposite high energy limit the summation can be replaced by integration over the continuous variable 
x = 27r/(eZ s )(n + 1/2), leading to 



G 



d(e) C s 



d N l s J xn l + x 3 el 



dx, Xq = — , (27) 



and finally 



d(e) C s (1 4 l_\ c hxC ± 
d N ir 62 1 \3> '3' zgj ~" 3^0 



^ = ^ £2 ^(± 1,1,-4) ->f^f as £^oo. (28) 



The values of the parameters of the length distribution function are given by ^ s = 1 and C s = 3\/3/(27r). This value 
of C is very close to C = 8/tt 2 of § giving a slope of 3\/3/8 ~ 0.6495 of d(e)/dN for small energies compared to 
2/tt« 0.6366 of §. 

In Fig. ^| we present a characteristic example of length distribution probabilities for a small channel width w <C a. 
A plateau in the length distribution function P(l) for I < 1 is clearly visible. Here the crossover from the asymptotic 
power law for long lengths to the plateau is slightly better modelled by P c (l) compared to P s (l)- 

We also calculated the density of states for individual square billiards starting from the semiclassical formula (||) 
without employing the assumption of an ergodic distribution of initial conditions on the SN-boundary. For individual 
square billiards this formula is more appropriate than ([?]) before averaging over Fermi energy since it does not involve 
the assumption of an ergodic distribution of initial conditions of trajectories on the SN-boundary. In square billiards 
trajectories with equal length between two hits with the channel lead are organised in families (for a discussion see 
p9|). Orbits with channel index a = 1, • • • , N hit the SN-boundary under an angle sin# a = air/(kpw). For each 
channel a there is a finite number of orbit families A a with different lengths L\ a . Each member of an orbit family 
carries a weight 6\ a with which it contributes to the density of states which is then given by 

m = d ~f^ ± £ S xMa ±s(^-(n+ M (29) 

a=l{\ a } n=0 \ ■* ' 

with X){a } 3\ a — w f° r the sum of weights of the members A Q of a single channel a. The lengths L\ a of contributing 
trajectory families and their weights <5A a can be very efficiently calculated by means of an algorithm involving a finite 
number of continued fractions as described in |2£j . Fig. || shows the result of a quantum mechanical calculation for 
a billiard with w/a — 1/3 and N = 12 channels. The solid line is a 20-point average over the quantum mechanical 
data. The dashed line represents the semiclassical result obtained from Eq. (p9|). It follows the general trend of 
the quantum mechanical result but is somewhat lower. Also plotted as dotted and dot-dashed lines are the formulas 
for the averaged quasiparticle excitation spectrum obtained from Eq. (|2^), ( pif ) and (|2^). The linear behaviour of 
the quasiparticle spectrum for small excitation energies based on the asymptotic length distribution function P(l) 
without plateau is plotted as the long dashed line. The average of the quantum mechanical result is in good agreement 
with the results of Eq. (^3|) , (|24|) and ( |26| ) . Both the numerical semiclassical result (dashed line) and the quantum 
mechanical result show however additionally a pronounced oscillation around the mean value d(e) = d^ at energies 
e > 1. We observed this trend also for other individual square billiards with different parameter values of w/a and 
N. The origin of these oscillations has not yet been identified and deserves further investigations. 

3.2 Rough billiards 

In this subsection we apply the semiclassical result (||) for the quasiparticle density of states to the square billiard 
with additional surface roughness. In the following we discuss the situation where the width w of the superconducting 
channel is much smaller than the length a of one billiard side. Since a trajectory randomises once it hits the rough 
billiard walls due to off-scattering in arbitrary direction with equal probability the motion becomes ergodic on a time 
scale much smaller than the mean time between two hits with the superconducting part of the boundary. The latter 
can also be viewed as the escape time of trajectories from the billiard if it was open along the channel lead. The 
esc ape^ pro bability P(L) of trajectories of length L from an open chaotic billiard is known to be given asymptotically 

P( L)= * xpf--M (30) 



with the above defined Thouless length Lt ||. It is related to the mean escape time by t csc = Lx/vp- Note that 
averaging over different realizations of the rough billiard is effectively taken into account by the introduction of 
the smooth length distribution function P(L) by which a probabilistic description on the classical level is achieved. 
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Numerical calculations confirm the form (||(]) of P(L) for long trajectories. Using ( |30| ) the average density of states 
can be calculated from Eq. (^|) by evaluating the delta function integral directly and summing over n. The result is 



2 cosh(f) , E hv F 

^VMr with e = , E T = — - 

sinh 2 (*) E T 2L T 



d(e)=d N {-) . 1 2 \l\ with e =—,E T = W7 L . (31) 



This result was also derived in || in the context of a chaotic billiard coupled to a superconducting lead of small width 
w. It remains valid for billiards with boundary roughness as long as the time scale t er g on which the classical motion 
in the billiard becomes ergodic is much smaller than the mean escape time r osc . 

The semiclassical prediction ( |3l| ) is compared with numerical quantum mechanical results in Fig. ^. A 20-point 
average was taken over the numerical results for different channel widths w (solid curve) . It is in very good agreement 
with the analytical formula ( |3l| ) (dotted curve) as well as an evaluation of Eq. (|^) with numerically calculated length 
distribution functions P(L) (dashed curve). The Thoulcss energy is inversely proportional to twice the escape time 
t csc of an electron trajectory since for a complete transversal of a path hitting the SN-boundary the electron part and 
the retracing as a hole trajectory have to be added for one full cycle. 

The exponential suppression of the spectral density for the rough billiard in contrast to the linear rise in energy for 
the square billiard without roughness has its origin in the asymptotic behaviour of the length distribution function 
P{L). The density of low- lying quasiparticle excitations is determined by very long orbits in order to fulfill the delta 
function condition in Eq. (|3|) in the limit E — > 0. As a consequence of the exponential tail of the length distribution 
( ^p| ) the density of states is exponentially suppressed above the Fermi energy. In contrast the algebraic tail of P(L) 
in the square billiard without roughness leads only to a linear suppression in energy. 

In a different approach the Hamiltonian of an (isolated) chaotic mesoscopic billiard was modelled by a GOE ensemble 
from random matrix theory and its coupling to the superconducting leads was then taken into account by means of a 
coupling matrix [^J^] . There it was shown that the quasiparticle density of states vanishes exactly below an energy of 
approximately E c < 0.6-EV (see dashed-dotted line in Fig. 4). As discussed at the end of Sec. II the difference between 
the two results can be traced back to the diagonal approximation employed in the semiclassical theory. Despite this 
approximation the quantum mechanical results of Fig. || are in good agreement with the semiclassical prediction of an 
exponential suppression of the density of states. This is due to the fact that the random matrix theory of H predicts 
a gap in the excitation spectrum in the limit of an infinity number of channels N — > oo (since it is essentially an 
expansion in the parameter 1 /N) . For a finite number of channels the behaviour of the spectral density around E c is 



expected to be smoothed out [Ell. 

3.3 Effect of a magnetic field in rough billiards 

In this subsection we consider the effect of a magnetic field on the density of quasiparticle states. The magnetic 
field B is uniform and points in the direction perpendicular to the billiard area. The superconducting lead itself is not 
penetrated by the flux. We study the perturbative regime of small magnetic fields where the cyclotron radius is much 
larger than the length scale a defined by the size of the mesoscopic billiard itself. The trajectories are then unchanged 
and the only effect of the magnetic field is an additional phase acquired by each trajectory which is proportional to 
the directed flux enclosed by the trajectory j3lj]. Formula (EH for the average density of states is then modified to 
include the flux dependent phase for each trajectory. The result is 

W 1 L(y,a) 

d(E)= d -fJdyJd(sina) J ds £ S (^M _ (n + ^ + 27r £|^) , (32 ) 



o 

where $o = ch/e is the flux quantum. Electron trajectories which traverse the billiard in opposite directions between 
two hits with the superconducting part of the boundaries acquire flux of same magnitude but opposite sign. This 
corresponds to a splitting of energy levels linear in B which are degenerate in the absence of a magnetic field |52| . 
For the generic case of rough or chaotic billiards a statistical description can be used in the same way as was done 
in the previous section for the magnetic field free case. In addition to the length distribution P(L), which remains 
unchanged, the distribution Pl(8) of the directed area O enclosed by the trajectories of given length L must be 
specified. For long trajectories it is given by a Gaussian of the form p0|.pl| 



PL(eH 7^ exp (~^)' 



(33) 
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where 8 is the directed area enclosed by the orbit. A heuristic argument for this distribution goes like follows: for 
a long orbit the length is on average proportional to the number of bounces with the rough billiard walls. At each 
bounce the trajectory is randomised and therefore the area swept between two successive bounces can be viewed 
as a random variable with zero mean value. Trajectory segments are independent of each other. The total area O 
accumulated by a trajectory is therefore a sum of independent random variables and its distribution is a Gaussian. It 
follows that the area distribution integrated over length has an exponential form: P(O) = (20) _1 cxp(— 0/0) with 
= ^L T a L /2. 

Using the above area distribution and an ergodic distribution of initial conditions along the superconducting channel 
boundary the density of states can be rewritten as 



d(E,B) = ^j- dLP(L)L d<3P L (Q) 

A JO J-oo 



J EL . U 2ttB 

x> 5[ (n+-)n-\ 

^ \hv F 1 2 7 $ 

Transforming the sum over delta functions by use of the Poisson formula and performing the integration over areas 
and lengths successively one finally gets the following expression for the average spectral density of quasiparticlc 
excitations as a function of energy and magnetic flux: 

dM-tiftg-ll' ^-W; , *=ff^/2^;. (34, 

The dimensionless quantity <fi is the effective flux through the billiard measured in units of the flux quantum $o = ch/e. 
This is evident when it is written as 4> = 47r_B0/$o- has the meaning of an effective area. The dimensionless flux 
can also be expressed as cj> = ^tot/^cr, where $ to t = BA is the total flux through the billiard and $ cr = A/ (47rO)<I>o. 
As will be discussed in the following in the regime <I> tot < $ cr ((f) < 1) the exponential suppression of the DOS 
near the Fermi energy persists while for larger fluxes the DOS converges towards a constant value. Fig. || shows 
the quasiparticle excitation spectrum as a function of energy for two different values of the flux parameter <fi. With 
growing flux the average density of states acquires the value d/v of the isolated billiard also for values E < Et- 



At the Fermi energy e = Eq. (34) can be summed and gives 

(35) 



sinh — 



7T , x f 
- coth - + 1 



Eq. ( |35| ) should not be understood in the sense of an expression for the density of states exactly at the Fermi level. 
Properly interpreted it is proportional to the quasiparticle density of states as a function of the magnetic flux when the 
density of states is averaged over an energy window < E < 5 between the Fermi energy and the mean level spacing 
S of the isolated billiard. (Note that the averaging is over an energy scale much smaller than Et-) It has been shown 
that even at fluxes ^ > 1 a minigap of order of the mean level spacing is present in the quasiparticle density of states 
at the Fermi level Therefore the quasiparticle DOS exactly at the Fermi level is always zero. According to 

Ref. [0 for a billiard whose isolated dynamics is described by a GUE matrix ensemble (which corresponds to (j> > 1) 
the DOS is given by gJgueC-E') / d>N = 1 — sin(27rc?Ari?)/(27rc?Ari?). When dau^iE) / 'dx is averaged over the interval 
< E < S it gives the value g — 0.77. Using this asymptotic value for large fluxes we therefore expect the DOS d$ 
averaged over the window < E < S to be given by 

d 5 (0» =gd{0,<t>). (36) 

The factor g has the effect that ds(Q, <j>) saturates at a value smaller than 1 in large magnetic fields, which is the 
consequence of the existence of the minigap. 

Fig. H compares the theory of Eq. ( |36"| ) (dashed line) with a numerical quantum mechanical calculation (solid line). 
It shows that ds(0, (f>) is exponentially suppressed on the scale (f> < 1. The quantum mechanical results were obtained 
for a rough billiard (sides a = 1 and b = 1.5) with 10 different values of the width between w — 0.2 and w = 0.38. In 
Ref. |3(| numerical evidence was put forward that the effective area scales like = ao^4 5 ^ 4 w -1 / 2 with the parameters 
of the billiard («o is a numerical parameter). The flux scale $ cr entering into the magnetic field dependence of the 
DOS then scales as $ cr = (4-Kao)~ 1 w 1 / 2 A~ 1 / 4 $>o- The quantum mechanical calculations for different channel widths 
w and fixed billiard area A confirm this scaling property. Determining ccq = 0.1 from the numerical data the flux 
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scale $ cr is fixed for each value of w. The quantum mechanical data points plotted as a function of the flux ^tot/^cr 
and its average (solid line) are in excellent agreement with the semiclassical theory of Eq. (^) . 

The difference between Eq. (|35|) which would lead to a finite density of states at e = at nonzero flux and the 
exact result of Ref. |i~7| which predicts the minigap demonstrates the importance of nondiagonal contributions in the 
trace formula (0) for energies E < S. As it is known the diagonal approximation, which we employed, breaks down on 
the scale of the mean level spacing [^6|,[53| and therefore semiclassics within the diagonal approximation does not lead 
to the correct result on the scale of the mean level spacing itself. When the density of states is averaged over a scale 
of the mean level spacing 6 or larger Eq. (|3^) describes the behaviour of the spectral density correctly as discussed 
above. 

Finally we would like to emphasize that the semiclassical theory of the proximity effect on the density of states in 
rough or chaotic Andreev billiards allows for a very intuitive picture of its destruction in the presence of a magnetic 
flux. The suppression of the density of states below Ex is a consequence of the coupling between electron and 
hole like quasi-particles at the SN-interface. Classically the return probability of trajectories to the SN-interface is 
given by P(L). Without magnetic field a path and its reverse interfere constructively in Eq. (S) for the density of 
states. In the presence of a magnetic flux the phase difference between a path and its reverse leads to destructive 
interference. Effectively the quantum mechanical return probability is reduced compared to its classical value P{L). 
Stated differently the coupling of the mesoscopic billiard to the superconducting lead (the possibility of electron hole 
conversion) is reduced in the presence of a magnetic flux and thus the proximity effect is suppressed. With growing 
magnetic flux the billiard will effectively look like an isolated mesoscopic conductor and its average density of states 
will therefore acquire the value d^. The effect is reminiscent of the weak localisation correction to the classical 
reflection coefficient |3(J and of enhanced orbital magnetism |3l]] in ballistic semiconductor microstructures. 



4 Summary 



We studied the quasiparticle excitation spectrum of a mesoscopic conductor modelled by a billiard in proximity 
to a superconducting lead. The expression for the density of states was derived from the semiclassical scattering 
matrix formulation. It was shown to be equivalent to the result derived previously from the Eilenberger equation for 
the quasiparticle Greens function when the diagonal approximation is applied to traces of powers of the scattering 
matrix. Applications to a square billiard as an example of a classical integrable system and a square billiard with 
surface roughness as an ergodic system were considered. For the square billiard without roughness it was shown that 
the deviation of the length distribution function of trajectories hitting the superconducting part of the boundary at 
small lengths is responsible for the saturation of the average density of states at large energies. For the rough billiard 
semiclassics predicts an exponential suppression of the density of states at energies smaller than the Thouless energy. 
The difference to the random matrix result of an energy gap of order of the Thouless energy can be traced back to the 
diagonal approximation employed in the semiclassical theory. The semiclassical predictions are in good agreement with 
quantum mechanical calculations for the integrable as well as for the ergodic billiard. Finally we derived expressions 
for the effect of a weak magnetic field on the density of states in the rough billiard. A magnetic field destroys the 
proximity effect on the density of states. We interpreted this fact as a phase phenomenon involving identical but 
reversed paths which hit the superconducting part of the billiard boundary. Due to destructive interference in the 
presence of a magnetic flux the quantum mechanical return probability of trajectories to the superconducting lead is 
smaller than its classical value effectively reducing the coupling between the billiard and the superconducting lead. 

We enjoyed discussions with A. Altland, C. Beenakker, K. Frahm, B. Mehlig and H. Schomerus and thank M. Sieber for a 
critical reading of the manuscript. WI would especially like to thank F. Mota-Furtado and P. F. O'Mahony for encouraging 
him to work on Andreev billiards and for discussions on the subject in the early stage of this work. 
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Fig. 1 
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FIG. 1. Square billiard with rough boundaries. Roughness is modelled by isotropic scattering when an electron 
or hole hits the normal billiard boundary (reflection in any direction with equal probability) in contrast to specular 
reflection at a smooth boundary. At the SN boundary electrons are retro-reflected (Andreev reflected) as holes and 
vice versa. 
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Fig. 2 
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FIG. 2. Histogram of the length distribution probability P(L/Lt) for a square billiard obtained by evaluating Eq. 
( p9| ) by the method of continued fractions. Parameters: Length of the billiard sides a — 1, width of the superconducting 
lead w — 0.1. The histogram is an average over different numbers of channels ranging from N = 96 to TV = 105. 
Dashed line: asymptotic power law for large I. Dotted line: P c (l), Eq. @. Dot-dashed line: P s (l), Eq. @. 
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Fig. 3 
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FIG. 3. Average density of states d(e)/djv as a function of the energy e = E/Et- The circles are quantum 
mechanical calculations for a square billiard of side length a = 75 and channel width w = 25 (N = 12). The solid 
line is a 20-point average over the set of data points. Dashed line: Semiclassical result from Eq. (^). Dotted line: 
DOS using P C {1) from Eq. (po[). Dot-dashed line: DOS using P 3 (l) from Eq. © Long dashed line: Low energy limit 
d(e)/d,N = (2/ir)e based on the probability distribution P(l) — C/l 3 . 
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Fig. 4 



1.5 




FIG. 4. Average density of states for a rough square billiard. Data points are quantum mechanical energy eigen- 
values for a billiard of side length a = 75 and channel widths w < 40. The solid line is a 20-point average of the 
numerical data. Dashed curve: semiclassical calculation based on Eq. (Q). Dotted curve: Analytical expression Eq. 
(30). Dot-dashed line: Random matrix result of Ref. [3]. 
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Fig. 5 
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Fig. 6 




FIG. 6. Spectral density ds(0,cj>)/dN of quasiparticle excitations as a function of the flux <f) = $tot/3>cr for the 
rough square billiard. The density of states is averaged over the small energy interval < E < S, where <5 is the mean 
level spacing of the isolated billiard. Data points: quantum mechanical calculation width 10 different channel width 
w. Solid line: Average over quantum mechanical data points. Dashed line: semiclassical theory, Eq. (B5I) and (p6|). 
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